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, Radial basis functions provide highly useful and flexible interpolants to mul- 

^— ( ' tivariate functions. Further, they are beginning to be used in the numerical 

(~| , solution of partial differential equations. Unfortunately, their construction 

requires the solution of a dense linear system. Therefore much attention has 
been given to iterative methods. In this paper, we present a highly efficient 
preconditioner for the conjugate gradient solution of the interpolation equa- 
tions generated by gridded data. Thus our method applies to the correspond- 
ing Toeplitz matrices. The number of iterations required to achieve a given 
\ tolerance is independent of the number of variables. 

CN ' 1. Introduction 

o 
o 



A radial basis function approximation has the form 

n 



^ ' where ip: [0, cxj) — >■ M is some given function, (yj)" are real coefficients, and 

^ . the centres {xj)^ are points in R'^; the norm || • || will be Euclidean throughout 

this study. For a wide class of functions ip, it is known that the interpolation 
matrix 

A = {ifiWxj - Xk\\))lk=i 

is invertible. This matrix is typically full, which fact has encouraged the 
study of iterative methods. For example, highly promising results have been 
published in the use of radial basis functions in collocation and Galerkin 
methods for the numerical solution of partial differential equations (see 
Franke and Schaback (1998) and Wendland (1999)), but direct solution 
limits their applicability to fairly small problems. The use of the precon- 
ditioned conjugate gradient algorithm was pioneered by Dyn, Levin and 
Rippa (1986), and some stunning results for scattered data were presented 
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recently in Faul and Powell (1999), although the rapid convergence described 
there is not fully understood. Therefore we study the highly structured case 
when the data form a finite regular grid. The conjugate gradient algorithm 
has been applied to Toeplitz matrices with some success; see, for instance, 
Chan and Strang (1989). However, since our matrices are usually not posi- 
tive definite and often possess elements that grow away from the diagonal, 
the preconditioners of Chan and Strang (1989) are not suitable. However, 
the matrices have the property that their inverses tractable more tractable. 
Specifically, the detailed study of the spectra of the associated Toeplitz op- 
erators presented in Baxter (1992) and Baxter (1994) allows us to create 
highly efficient preconditioners by inverting relatively small finite sections 
of the bi-infinite symmetric Toeplitz operator, and this construct is also 
easily understood via Toeplitz theory. 

Let n be a positive integer and let An be the symmetric Toeplitz matrix 
given by 

An = {^U-k))l,^_^, (1.1) 

where (/9:M — t- M is either a Gaussian (fix) = exp(— Ax^) for some positive 
constant A) or a multiquadric {(p{x) = {x'^ + c^)^/'^ for some real constant 
c). In this paper we construct efficient preconditioners for the conjugate 
gradient solution of the linear system 

Anx = f, /eM2"+\ (1.2) 

when (/9 is a Gaussian, or the augmented linear system 

AnX + ey = f, (1.3) 
e^x = 0, (1.4) 

when ip is a multiquadric. Here e = [1, 1, . . . , 1]"^ € R^"^-'^ and y G M. Section 
2 describes the construction for the Gaussian and Section 3 deals with the 
multiquadric. Of course, we exploit the Toeplitz structure of An to perform 
a matrix- vector multiplication in O(nlogn) operations whilst storing 0{n) 
real numbers. Further, we shall see numerically that the number of iterations 
required to achieve a solution of (|1.2p or (|1.4p to within a given tolerance 
is independent of n. The Matlab software used can be obtained from my 
homepage. 

Our method applies to many other radial basis functions, such as the in- 
verse multiquadric {(p{x) = (x^ -|- c^)"-"^/^) and the thin plate spline {(p{x) = 
x^log|x|). However, we concentrate on the Gaussian and the multiquadric 
because they exhibit most of the important features of our approach in a 
concrete setting. Similarly we treat the one-dimensional problem merely to 
avoid complication; the multidimensional case is a rather slight generaliza- 
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tion of this work. Let us remark that the analogue of (jl.ip is the operator 

4'^) = ((^(j-A;))^.,g[_„^„],, (1.5) 

and we shall still call a Toeplitz matrix. Moreover the matrix-vector 
multiplication 

4')^= I E M-k\\)xk] , (1.6) 

where || • || is the Euclidean norm and x = (xj)jg[_„ can still be cal- 
culated in 0{N log N) operations, where N = {2n + l)'^, requiring 0{N) 
real numbers to be stored. This trick is a simple extension of the Toeplitz 
matrix- vector multiplication method when d = 1. 

2. The Gaussian 

It is well-known that the Gaussian generates a positive definite interpola- 
tion matrix, and its functional decay is so rapid that preconditioning the 
conjugate gradient algorithm is not necessary. However, it provides a use- 
ful model problem that we shall describe here before developing the ideas 
further in the following section. 

Our treatment of the preconditioned conjugate gradient (PCG) method 
follows Section 10.3 of Golub and Van Loan (1989), and we begin with 
a general description. We let n be a positive integer and A G R"^" be an 
arbitrary symmetric positive definite matrix. For any nonsingular symmetric 
matrix P € M"^" and b E we can use the following iteration to solve the 
linear system PAPx = Ph. 

Algorithm 2.1. Choose any xq in M". Set tq = Pb — PAPxq and do = rg. 
For = 0, 1, 2, . . . do begin 
Ofc = rlrk/dlPAPdk 
Xk+i = Xk + akdk 
rk+i = rfc - QkPAPdk 
bk = rl^iTk+i/rln 

4+1 = rk+i + bkdk 
Stop if or 114+1 II is sufficiently small, 

end. 

hi order to simplify Algorithm 12. II define 

C = P"^, Ck = Pxk, Tk = Ppk and 5k = Pd^. (2.1) 
Substituting in Algorithm 12.11 we obtain the following method. 
Algorithm 2.2. Choose any in l^"- Set po = b — A^o, 6q = Cpo- 
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For k = 0,1,2, .. . do begin 
«fe = p1Cpk/5lA5k 

ik+l = ik + O-k^k 

Pk+i = Pk- akA5k 
h = pl+iCpk+i/plCpk 

^k+l = Cpk+i + bkSk 
Stop if llpfe+ill or is sufficiently small, 

end. 

It is Algorithm 12.21 that we shall consider as our PCG method in this 
section, and we shall call C the preconditioner. We see that the only re- 
striction on C is that it must be a symmetric positive definite matrix, but 
we observe that the spectrum of CA should consist of a small number of 
clusters, preferably one cluster concentrated at one. At this point, we also 
mention that the condition number of C A is not a reliable guide to the ef- 
ficacy of our preconditioner. For example, consider the two cases when (i) 
CA has only two different eigenvalues, say 1 and 100, 000, and (ii) when CA 
has eigenvalues uniformly distributed in the interval [1,100]. The former 
has the larger condition number but, in exact arithmetic, the answer will be 
achieved in two steps, whereas the number of steps can be as high as n in 
the latter case. Thus the term "preconditioner" is sometimes inappropriate, 
although its usage has become standard. 

In this paper we concentrate on preconditioners for the Toeplitz matrices 
generated by radial basis function interpolation on a (finite) regular grid. 
Accordingly, we let A be the matrix An of (jl.ip and let (p{x) = exp(— a;^). 
Thus An is positive definite and can be embedded in the bi-infinite symmet- 
ric Toeplitz matrix 

^oo = {^{j - A:)),, fce^ . (2.2) 

The classical theory of Toeplitz operators (see, for instance, Grenander and 
Szego (1984)) and the work of Baxter (1994) provide the relations 

Sp An C Sp Aoo = k(7r), a(0)] C (0, oo), (2.3) 

where a is the symbol function 

^^(0 = j;<^(e + 2^/c), CeR, (2.4) 
kez 

and Sp denotes the spectrum of the operator A^o- Further, Theorem 9 
of Buhmann and Micchelli (1991) allows us to conclude that, for any fixed 
integers j and k, we have 

hm {A-')j^k = {A^\k. (2.5) 

n— >-oo 

It was equations (|2.3p and (12. 5j) which led us to investigate the possibility 
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of using some of the elements of for a relatively small value of n to 
construct preconditioners for A]sf, where may be much larger than n. 
Specifically, let us choose integers < m < n and define the sequence 

= (^n^)jo, j = -m, . . . , m. (2.6) 



We now let Cat be the (2A^ + 1) x (2A^ + 1) banded symmetric Toeplitz 
matrix 



/ Co 



Cat 



V 



(2.7) 



Co 



We claim that, for sufficiently large m and n. Cat provides an excellent 
preconditioner when A = Ai\[ in Algorithm 12. 2[ Before discussing any theo- 
retical motivation for this choice of preconditioner, we present an example. 
We let n = 64, m = 9 and N = 32, 768. Constructing An and calculating 
the elements {{A~^)jQ : j = 0, 1, . . . , m} we find that 

^ 1.4301 X 10° \ 
-5.9563 X 10-1 
2.2265 X 10-^ 
/co\ -8.2083 x 10-2 
Cl 3.0205 X 10-2 

-1.1112 X 10-2 
\cgj 4.0880 x 10-3 

-1.5039 X 10-3 
5.5325 X 10-"^ 
V -2.0353 X 10- V 



(2. 



by 



Now Cn can be embedded in the bi-infinite Toeplitz matrix Cqo defined 



Cj-k, \j -k\ <m, 
0, \j - k\> m, 



{2S 



{Coo)jk " 

and the symbol for this operator is the trigonometric polynomial 

m 

^c^iO= E ^e^- (2.10) 



j=~rn 



In Figure [2?T] we display a graph of (Tc^ for < ^ < 2tt, and it is clearly 
a positive function. Thus the relations 



Sp Cn C Sp Coo = WcAO : e G [0, 27r]} C (0, oo) 



(2.11) 
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1. 

1.6- 

1.4- 

1.2- 

1 - 

0.8 - 

0.6 

0.4- 

0.2 - 




Fig. 2.1. The symbol function for Co, 



imply that Cn is positive definite. Hence it is suitable to use Cat as the 
preconditioner in Algorithm 12.21 Our aim in this example is to compare this 
choice of preconditioner with the use of the identity matrix as the precondi- 
tioner. To this end, we let the elements of the vector h of Algorithm 12.21 be 
random real numbers uniformly distributed in the interval [—1, 1]. Applying 
Algorithm 12.21 using the identity matrix as the preconditioner provides the 
results of Table 12.11 Table 12.21 contains the analogous results using (12. 7p and 
p.Sp . In both cases the iterations were stopped when the residual vector 
satisfied the bound ||rfe+i ||/||6|| < 10^^^. The behaviour shown in the tables 
is typical; we find that the number of steps required is independent of N 
and b. 

Why should (|2.7p and (12. Sp provide a good preconditioner? Let us consider 
the bi-infinite Toeplitz matrix Coo^oo- The spectrum of this operator is 
given by 



SpC 



A 

oo^oo 



{^C^(eMO:CG[0,2^]}, 



(2.12) 



where a is given by (12. 4p and ctc^ by (I2.10p . Therefore in order to concen- 
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Table 2.1. No preconditioning 



Iteration 


Error 


1 


2.797904 X 10^ 


10 


1.214777 X 10-2 


20 


1.886333 X 10-^ 


30 


2.945903 X 10"^" 


33 


2.144110 X 10"" 


34 


8.935534 x lO^^^ 



Table 2.2. Using (|2.7p and (j2.8p as the preconditioner 



Iteration 


Error 


1 


2.315776 X 10-1 


2 


1.915017 X 10"^ 


3 


1.514617 X lO"'^ 


4 


1.365228 X 10-" 


5 


1.716123 X 10-^"^ 



trate Sp Coo^oo at unity we must have 

^C^(eMO«l, CG[0,27r]. (2.13) 

In other words, we want crc^ to be a trigonometric polynomial approximat- 
ing the continuous function l/cr. Now if the Fourier series of 1/cr is given 

by 

^~HO = T.lje''^, eeK, (2.14) 

then its Fourier coefficients {'~fj)j^z are the coefficients of the cardinal func- 
tion X for the integer grid, that is 

xix) = Y.lM^-j)^ XGM, (2.15) 

and 

X{k) = Sok, keZ. (2.16) 

(See, for instance, Buhmann (1990).) Recalling (|2.5p . we deduce that one 
way to calculate approximate values of the coefficients {^j)jez is to solve 
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the hnear system 



where e° = (5io)"=-n G M^^+i. We now set 



< j < m, 



(2.17) 



(2.18) 



and we observe that the symbol function a for the Gaussian is a theta func- 
tion (see Baxter (1994), Section 2). Thus cr is a positive continuous function 
whose Fourier series is absolutely convergent. Hence 1/cr is a positive con- 
tinuous function and Wiener's lemma (Rudin (1973)) implies the absolute 
convergence, and therefore the uniform convergence, of its Fourier series. 
We deduce that the symbol function can be chosen to approximate 1/a 
to within any required accuracy. More formally we have the 

Lemma 2.1. Given any e > 0, there are positive integers m and no such 
that 



1 



< e, 



[0,27r], 



(2.19) 



j=-m 



for every n> uq, where c^") = (c^ )"^_^ is given by (j2.17p . 

Proof. The uniform convergence of the Fourier series for implies that 
we can choose m such that 



^(0 J2 ^^^''^ - 1 



< e, 



CG [0,27r]. 



j=-m 



By (|2.5p . we can also choose no such that max{|7j — c! 
e, when n > hq. Then we have 



-m, 



(2.20) 

, m} < 



^(0 E 



j=-m 



j=—m j=—m 

< e[l + (2m+l)||(T||oo]. 

Since e is arbitrary the proof is complete. 

3. The Multiquadric 

The multiquadric interpolation matrix 

A = (^{\\xj -Xk\\)) _ , 



(2.21) 

□ 



(3.1) 
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where ip{r) = (r^ + c^)-^/^ and (xj)^^-^ are points in R'^, is not positive 
definite. In Micchelli (1986), it was shown to be almost negative definite, 
that is for any real numbers {yj)^^i satisfying ^ y j = we have 

n 

X] VjVk^iW^j - ^fcll) ^ 0- (3-2) 
j,k=i 

Furthermore, inequality (13. 2p is strict when n > 2, the points (xj)"^^ are 
all different, and ^ \yj\ > 0. In other words, A is negative definite on the 
subspace (e)-*-, where e = [1, 1, . . . , 1]^ G M". 

Of course we cannot apply Algorithms 12.11 and 12.21 in this case. However, 
we can use the almost negative definiteness of A to solve a closely related 
linearly constrained quadratic programming problem: 

minimize —S/^A^ — b^^ 
subject to e^S, = 0, 

(3.3) 

where b can be any element of M". Standard theory of Lagrange multipliers 
guarantees the existence of a unique pair of vectors ^* € M" and r]* £ 
satisfying the equations 

AC + er]* = b, 

e^C = 0, (3.4) 

where rj* is the Lagrange multiplier vector for the constrained optimization 
problem (13. 3p . We do not go into further detail on this point because the 
nonsingularity of the matrix 

is well-known (see, for instance, Powell (1990)). Instead we observe that one 
way to solve p.4p is to apply the following modification of Algorithm 12. II to 
(133]). 

Algorithm 3.1. Let P be any symmetric nx n matrix such that kerP = 
(e). 

Set xq = 0, ro = Pb — PAPxq, cIq = tq. 
For /c = 0, 1, 2, . . . do begin 
Ofc = rlrk/dlPAPdk 
Xk+i = Xk+ akdk 
Tfc+i = rfc - akPAPdk 
bk = rl_^irk+i/rlrk 
dk+i = rk+i + bkdk 
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Stop if ||r(fc+i|| or ||(i/!c+i|| is sufficiently small, 
end. 

We observe that Algorithm lS.ll solves the linearly constrained optimization 
problem 

minimize —x^PAPx — b^Px 
2 



subject to e^x = 0. 



(3.6) 



Moreover, the following elementary lemma implies that the solutions ^*of 
(13. 4p and x* of (13. 6|) are related by the equations ^* = Px*. 

Lemma 3.1. Let S be any symmetric n x n matrix and let K = kerS". 
Then S : is a bijection. In other words, given any h G there 

is precisely one a G such that 

Sa = b. (3.7) 



Proof. For any n x n matrix M we have the equation 

R" = kerMelm M^. 
Consequently the symmetric matrix S satisfies 

= ker 5 e Im S, 

whence Im S = K^. Hence for every b € there exists a € M" such 
that Sa = b. Now we can write a = a + /3, where a € and /3 ^ K are 
uniquely determined by a. Thus Sa = Sa = 6, and (|3.7p has a solution. If 
a' G also satifies ()3.7p , then their difference a — a' lies in the intersection 
K n K-^ = {0}, which settles the uniqueness of a. □ 

Setting P = S and K = (e) in Lemma 13.11 we deduce that there is exactly 
one X* G (e)^ such that 

PAPx* = Pb, 

and PAP is negative definite when restricted to the subspace (e)^. Follow- 
ing the development of Section 2, we define 

C = P\ (k = Pxk, and dk = Pdk, (3.8) 

as in equation (12. ip . However, we cannot define pk by (12. ID because P is 
singular. One solution, advocated by Dyn, Levin and Rippa (1986), is to 
use the recurrence for (pk) embodied in Algorithm 12.11 without further ado. 

Algorithm 3.2. Choose any in (c)^- Set pQ = b — A^q and 5q = Cpo- 



Preconditioned CG and RBFs 



11 



For k = 0,1,2, .. . do begin 
Ofc = plCpk/djASk 
Cfc+i = Cfc + ak^k 
Pk+i = Pk — akA6k 
h = pl+^Cpk+i/plCpk 
^k+i = Cpk+i + hkSk 
Stop if ll/Ofe+ill or is sufficiently small, 

end. 

However this algorithm is unstable in finite precision arithmetic, as we 
shall see in our main example below. One modification that successfully 
avoids instability is to force the condition 

Pk e (e)^ (3.9) 

to hold for all k. Now Lemma [3.11 implies the existence of exactly one vector 
Pk G (e)"*" for which Ppk = Vk- Therefore, defining Q to be the orthogonal 
projection onto (e)-*-, that \s Q : x ^ x — e{e^ x) / {e^ e) , we obtain 

Algorithm 3.3. Choose any in (e)^- Set po = — ^Co); ^^o = CpQ. 
For /c = 0, 1, 2, . . . do begin 

a-k = plCpk/5lA5k 

ik+l = ik+ CLkSk 

Pk+i = Q{pk - akA6k) 
h = pl^iCpk+i/plCpk 
^k+i = Cpk+i + bk5k 
Stop if i|pfc+i|| or pfc+ill is sufficiently small, 
end. 

We see that the only restriction on C is that it must be a non-negative 
definite symmetric matrix such that kerC = (e). It is easy to construct 
such a matrix given a positive definite symmetric matrix D by a rank one 
modification: 

The Cauchy-Schwarz inequality implies that x'^Cx > with equality if 
and only if x G (e). Of course we do not need to form C explicitly, since 
C : X I— 7> Dx — {e^ Dx/e^ De)De. Before constructing D we consider the 
spectral properties of A^o = {^{j — k))j^kez in more detail. 

A minor modification to Proposition 5.2.2 of Baxter (1992) yields the 
following useful result. Let us say that a complex sequence {yj)z is zero- 
summing if it is finitely supported and satisfies Yl Vj — 0- The symbol 
function 

f^(e) = ^^^^(^ + 2^/^), CGM, (3.11) 
kez 
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now requires the distributional Fourier transform of the multiquadric. In 
the univariate case, this is given by 

(^(0 = -(2c/|e|)i^i(c|e|), CGM\{0}, (3.12) 

where Ki is a modified Bessel function. The symbol function is studied 
extensively in Baxter (1994). 



25 



20- 



15 - 



10 



5 - 








1 2 3 4 5 6 7 



Fig. 3.2. The reciprocal symbol function 1/a for the multiquadric. 



Proposition 3.2. For every r] G (0, 2tt) we can find a set {(yj )jez '■ n = 
1, 2, . . .} of zero-summing sequences such that 

i™, E - \yt^\" = (3-13) 

Proof. We adopt the proof technique of Proposition 5.2.2 of Baxter (1992). 
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For each positive integer n we define the trigonometric polynomial 

n-l 
fc=0 

and we recall from Section 2 of Baxter (1994) that 

Kn{i)='^^^:^ = \Lnm\ (3.14) 

n sm 4/2 

where Kn is the nth degree Fejer kernel. We now choose {y^^^)j£z to be the 
Fourier coefficients of the trigonometric polynomial Ln{i — r?)sin^/2, 
which implies the relation 

|^yf)e*^-«|' = sin2^/2if„(e-??), 

and we see that {y^-^^)jez is a zero-summing sequence. By the Parseval 
relation we have 

, , r2TT 

^|yf)|2 = (2vr)-i / sin2e/2i^„(e-7?)de (3.15) 

and the approximate identity property of the Fejer kernel (Zygmund (1988), 
p. 86) implies that 

f27r 



sin^ r//2 = lim (27r)"^ / sin^ C/2 K^iC - r]) = lim V \y\ 



{")|2 
J 

iG2 



Further, because a is continuous on (0,27r) (Baxter (1994), Section 4.4), we 
have 

r-2n 

sm^7]/2 a{r]) = lim (27r)"W sin^ S,/2 Kn{S, - ri)a{C) d^ 

n-+oo _/q 



lim V vy'vl^^'ifij - k). 



□ 



Thus we have shown that, just as in the classical theory of Toeplitz op- 
erators (Grenander and Szego (1984)), everything depends on the range of 
values of the symbol function a. Because a inherits the double pole that ip 
enjoys at zero, we have a: (0, 27r) i-> (cr(7r), oo). In Figure [32] we display the 
function (t~^ . 

Now let m be a positive integer and let {dj)JL_^ be an even sequence 
of real numbers. We define a bi-infinite banded symmetric Toeplitz matrix 
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-Coo by the equations 

Thus (-Doo-4oo)jfc = ipij - k) where tp{x) = YdL-m^^i^i^ ~ 0- Further 

E ymij-k) = {2Ti)-^ r\Y.y^^''^f^D^(^)^{Odc, (3.17) 

j,k€Z jeZ 

where the symbol function croac fo^" the Toephtz operator Doo is given by 

m 

^D^(0= E rf.e^'^eeK. (3.18) 

j=-m 

Now the function aaD^ is continuous for G (0, 27r), so the argument of 
Proposition 13.21 also shows that, for every rj G (0, 27r), we can find a set 
{{y\"^^)jez : n = 1,2, . . . } of zero-summing sequences such that 



' fa^^ = (TD^{ri)a{rii). 3.19 

A good preconditioner must ensure that {cDoo(0^(0 • ^ ^ (Oj^tt)} is a 
bounded set. Because of the form of ctDoo have the equation 

E (ij- = 0. (3.20) 

3=-m, 

Moreover, as in Section 2, we want the approximation 

^z?oo(e)c^(0 eG(0,27r), (3.21) 

and we need ctDoo be a non-negative trigonometric polynomial which is 
positive almost everywhere, which ensures that every one of its principal 
minors is positive definite. 
Let us define 



(n) 



-(^n'),o' J = -m,...,m, (3.22) 
and 

<y-\i) = Y.^je''^, eeM. (3.23) 

Then Theorem 9 of Buhmann and Micchelli (1991) states that 

lim4"^=7j, (3.24) 
for any given fixed integer j. We shall use this fact to construct a suitable 
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Table 3.3. Preconditioned CG ~ m = 9, n = 64, TV = 2, 048 



Iteration I Error 



1 


3.975553 x 10^ 


2 


8.703344 X 10 


-1 


3 


2.463390 X 10 


-2 


4 


8.741920 X 10 


-3 


5 


3.650521 X 10 


-4 


6 


5.029770 X 10 


-6 


7 


1.204610 X 10 


-5 


8 


1.141872 X 10 


-7 


9 


1.872273 X 10 


-9 


10 


1.197310 X 10 


-9 


11 


3.103685 X IQ- 


11 



Table 3.4. Preconditioned CG - m = 9, n = 64, TV = 32, 768 



Iteration Error 



1 


2.103778 X 10^ 


2 


4.287497 x 10° 


3 


5.163441 X 10 


-1 


4 


1.010665 X 10 


-1 


5 


1.845113 X 10 


~3 


6 


3.404016 X 10 


-3 


7 


3.341912 X 10 


-5 


8 


6.523212 X 10 


-7 


9 


1.677274 X 10 


-5 


10 


1.035225 X 10 


-8 


11 


1.900395 X 10- 


10 



aooo- First we subtract a multiple of the vector [1, . . . , 1]'^ G M^™"''^ from 
(c^"^)J^_„ to form a new vector satisfying ^ d j = 0, and we 

observe that, by (j3.24p . cr£)^(^) is one-signed for all sufficiently large values 
of n. For the numerical experiments here, we have chosen n = 64 and m = 9. 
Thus, given 

AN=i^{j-k)] 

\ / j^k=-N 

for any N > n, we let Djy be any (2N + 1) x (2A^ + 1) principal submatrix 
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Table 3.5. Preconditioned CG ~ m = 1, n = 64, TV = 8, 192 



Iteration I Error 



1 


2.645008 X 10^ 


10 


8.632419 X 10° 


20 


9.210298 X 10 


-1 


30 


7.695337 x 10 


-1 


40 


3.187051 X 10 


-5 


50 


5.061053 X 10 


-7 


60 


7.596739 x 10 


-9 


70 


1.200700 x 10- 


10 


73 


3.539988 x 10" 


11 


74 


1.992376 X 10- 


11 



of Doo and define the preconditioner Cat by the equation 

= (3.25) 

D]\ie 

where e = [1, . . . , 1]^ € M^^*^^. We reiterate that we actuahy compute the 
matrix-vector product Cnx by the operations x i-> D]\!X—{e^ Dj^x/e^ DNe)e 
rather than by storing the elements of Cat in memory. 

Cat provides an excellent preconditioner. Tables 13.31 and 13.41 illustrate its 
use when Algorithm 13.31 is applied to the linear system 

Anx + ey = b, 
e^x = 0, 

(3.26) 

when iV = 2, 048 and = 32, 768 respectively. Here y G M, e = [1, . . . , 1]^ G 
]^27V+i g^^^ ^ g ]^2Af+i consists of pseudo-random real numbers uniformly 
distributed in the interval [—1,1]. Again, this behaviour is typical and all 
our numerical experiments indicate that the number of steps is independent 
of N. We remind the reader that the error shown is but that the 

iterations are stopped when either ||pfc+i!| or ||5fc+i|| is less than lO^^'^H^H, 
where we are using the notation of Algorithm 13.31 

It is interesting to compare Table [331 with Table [331 Here we have chosen 
m = 1, and the preconditioner is essentially a multiple of the second divided 
difference preconditioner advocated by Dyn, Levin and Rippa (1986). In- 
deed, we find that do = 7.8538 and di = d-i = —3.9269. We see that 
its behaviour is clearly inferior to the preconditioner generated by choosing 
m = 9. Furthermore, this is to be expected, because we are choosing a 
smaller finite section to approximate the reciprocal of the symbol function. 
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Fig. 3.3. The spectrum of CnAn for m = 1 and n = 64. 



However, because is a multiple of sin^^/2, this preconditioner still 

possesses the property that {crDadO'^iO ■ C ^ (Oj^tt)} is a bounded set of 
real numbers. 

It is also interesting to compare the spectra of CnAn for n = 64 and 
m = 1 and m = 9. Accordingly, Figures 13.31 and 13.41 display all but the 
largest nonzero eigenvalues of CnAn for m = 1 and m = 6 respectively. The 
largest eigenvalues are 502.6097 and 288.1872, respectively, and these were 
omitted from the plots in order to reveal detail at smaller scales. We see 
that the clustering of the spectrum when m = 9 is excellent. 

The final topic in this section demonstrates the instability of Algorithm 
13.21 when compared with Algorithm 13.31 We refer the reader to Table 13. 6^ 
where we have chosen m = 9, n = N = 64, and setting b = [1, 4, 9, ... , N'^]'^. 
The iterations for Algorithm 13. 3( displayed in Table 13.61 were stopped at 
iteration 108. For Algorithm [321 iterations were stopped when either || 
or became smaller than 10~^^||5||. It is useful to display the norm 

of 1 1 5^11 rather than ||pfe|| in this case. We see that the two algorithms 
almost agree on the early interations, but that Algorithm 13.21 soon begins 
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Fig. 3.4. The spectrum of C„^„ for m = 9 and n = 64. 



cycling, and no convergence seems to occur. Thus when can leave the 
required subspace due to finite precision arithmetic, it is possible to attain 
non-descent directions. 
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Table 3.6. Algorithms 3.3a & b - m = 1, n = 64, = 64, 
6=[l,4,...,iVr. 



Iteration 


ll^fcll - 3.3a 


ll^ll -3.3b 


1 


4.436896 x 10 


4.436896 x 10 


2 


2.083079 X 10^ 


2.083079 x 102 


3 


2.339595 x 10° 


2.339595 x 10° 


4 
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1.206041 x 10-1 


5 
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6 
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7 
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9.254943 x IQ-'^ 


8 
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9 


3.453789 x 10-^ 




10 


1.914126 X 10-3 




20 


4.628447 x 10-^ 




30 


3.696474 x 10-° 




40 


8.061922 x 10+3 




50 


2.155310 X 10° 




100 


3.374467 x 10-^ 
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